function intensity_ch_norm = intensity_channel_sim(img, pars)

%% input parameters
img_size_x = pars.img_size_x;
img_size_y = pars.img_size_y;
Nc = pars.Nc;
Ns = pars.Ns;
delta = pars.delta;
c = pars.c;
pyr_filt_coeff = pars.pyr_filt_coeff;
resample_filt_coeff = pars.resample_filt_coeff;
levels = pars.levels;
out_scale = pars.outscale;
dec_factor = pars.dec_factor;

output_nr = img_size_y/(dec_factor^out_scale);  %size of output conspicuity map
output_nc = img_size_x/(dec_factor^out_scale);

cs_ch = zeros(Nc*Ns,2);
intensity_ch = zeros(output_nr,output_nc);

% the center-surround levels for each channel
for i = 1:Nc
    clev = c(i);
    for j = 1:Ns
        slev = clev + delta(j);
        k = (i-1)*Ns+j;
        cs_ch(k,1) = clev;
        cs_ch(k,2) = slev;
    end
end
        
%% computation begins
pyr_select = 0; % 0 : gaussian pyramid, 1 : laplacian pyramid
pyr = pyramid_sim(img, pyr_filt_coeff, levels, pyr_select);         %Gaussian Pyramid generation

for lpi = 1:Nc*Ns   
    
    clev = cs_ch(lpi,1); %center channel for csd
    slev = cs_ch(lpi,2); %surround channel for csd
    cs_factor = 2^(slev-clev);
    
    center = pyr{clev+1};   %+1 to for matlab array index
    surround = pyr{slev+1};
    
    csd_out = csd_sim(center, surround, cs_factor);     %Center-surround difference

%     norm_out = maxnorm_sim(csd_out);                    % Max normalizer
    [norm_out, local_max_avg, norm_factor, lm_count, global_max] = maxnorm_sim_2(csd_out);
    
    resample_lev = abs(out_scale-clev);  % num of levels to output scale
    resample_out = resample_sim(norm_out, resample_filt_coeff, resample_lev); %Resample to output scale
    
    intensity_ch = intensity_ch + resample_out;    %Across-scale addition
end

% intensity_ch_norm = maxnorm_sim(intensity_ch);      %Final normalization
intensity_ch_norm = intensity_ch;      %Final normalization

